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Abstract 

The two-dimensional Ising model with nearest-neighbor ferromagnetic and long-range dipolar 
interactions exhibits a rich phase diagram. The presence of the dipolar interaction changes the 
ferromagnetic ground state expected for the pure Ising model to a series of striped phases as a func- 
tion of the interaction strengths. Monte Carlo simulations and histogram reweighting techniques 
applied to multiple histograms are performed to identify the critical temperatures for the phase 
transitions taking place for stripes of width /i = 2 on square lattices. In particular, we aim to study 
the intermediate nematic phase, which is observed for large lattice sizes only. For these lattice sizes, 
we calculate the critical temperatures for the striped-nematic and nematic-tetragonal transitions, 
critical exponents, and the bulk free-energy barrier associated with the coexisting phases. We 
also evaluate the long-term correlations in our time series near the finite-size critical points by 
studying the integrated autocorrelation time r as a function of the lattice size. This allows us to 
infer how severe the critical slowing down for this system with long-range interaction and nearby 
thermodynamic phase transitions is. 
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I. INTRODUCTION 



Simple two-dimensional (2D) Ising-like models may present a complex behavior with 
rich phase diagrams. A complex behavior may appear even in face of simple competing 
interactions. In general, a complex behavior results from the frustration phenomenon caused 
by these competing interactions at different length scales. 

The 2D Ising model with nearest- neighbor ferromagnetic exchange and dipolar interaction 
has been the focus of considerable interest because of its magnetic properties [1]. The 
Hamiltonian model is defined by the Ising lattice model with nearest-neighbor exchange 
interaction J > and dipolar contribution of strength g > 0, which is rewritten as 

<i,j> i<j «i 

where 5 = J/g. Here, we have adopted the convention [2j of summing over all distinct pairs 
of lattice spins at distances to define the constant g. The distances rjj are measured in 
units of lattice. 

The competition between the exchange and dipolar interactions results in a phase diagram 
that is rather sensitive to the 6 = J/g ratio [3l H]. For 6 < 0.4403, the model presents anti- 
ferromagnetic ground-state characterized by stable checkerboard like configurations [4]. For 
larger S values, the ground state changes to striped phases. These phases are characterized 
by magnetic domains displayed in stripes of alternating spins perpendicular to the lattice 
plane. Each stripe corresponds to a pattern of upward or downward spins, whose width h 
increases with 5 [Sj E]. As the temperature increases, for fixed couplings, an order-disorder 
transition takes place until the system reaches its usual paramagnetic phase. 

This theoretical aspect has been observed in experiments with ultrathin (i.e., a few atomic 
layers) metal films on metal substrates [H [7] below certain temperatures as a consequence 
of a reorientation transition of their spins. Apart from short-range interactions, long-range 
dipolar interactions are also present between the magnetic moments in these systems. As a 
matter of fact, even for the simple ferromagnetic Ising model, one should always include a 
long-range repulsive interaction between the magnetic moments in addition to the short- 
range exchange interaction term. Theoretical studies of magnetic systems have avoided the 
inclusion of the relatively small dipolar interaction compared with the exchange interaction. 
However, as we have mentioned the dipolar interaction plays an important role for quasi 
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2D-magnetic systems in characterizing metal films [T]. Furthermore, from a theoretical 
viewpoint, a long-range repulsive interaction for d-dimensional systems, with d < s < 
d + 1, in addition to the usual nearest-neighbor ferromagnetic interaction, leads to the 
absence of ferromagnetism for all temperatures [6l [8] . 

Efforts have been made to determine the phase diagram as a function of the coupling S 
from mean field approximation and numerical simulations [21 [3l El [9] . Analytical and Monte 
Carlo (MC) calculations have shown the following picture: For a fixed coupling 6 > 0.4403, 
the phase diagram initially exhibits the alternating striped spin configurations as described 
above, characterizing the so called smectic striped phase. As the temperature increases, the 
initial picture describes a transition into the tetragonal phase [ini [HI [El [13] characterized 
by states with orient at ionally disordered stripes but still preserving some of this structural 
form, which tends to the completely disordered paramagnetic state. Recent numerical results 
related to this striped-tetragonal phase transition indicate a continuous transition for /i = 1 
[H], a clear first-order transition for /i = 2 [TTl HI], a likely weaker first-order transition for 
h = 3 [H], and a continuous transition for h = 4 [2] and h = 8 [12]. On the other hand, 
a theoretical approach developed by Abanov et al. [I5] , based on a continuum limit for the 
Hamiltonian, has predicted a new domain phase in the phase diagram. In addition to the 
striped and tetragonal phases, there would be an intermediate Ising nematic phase, whose 
existence would depend on the parameters of the model. 

Different methods have been used to study the criticality of this model. As summarized 
above, the order of the phase transition between the low temperature smectic striped and the 
high temperature tetragonal phases is still under debate and seems to depend on the coupling 
6. Simulations have been hampered because large lattice size simulations are very CPU 
timeconsuming due to the dipolar term. In spite of these limitations, recent literature works 
have presented very satisfactory results, including numerical evidence about the intermediate 
Ising nematic phase [21 [TB]. The conclusion seems to be that the nematic phase can be 
observed in a narrow range of temperatures and is 5-dependent. This amounts to observing 
not only a single peak in the specific heat, initially corresponding to the expected striped- 
tetragonal phase transition, but a richer phase diagram comprising the striped-nematic and 
the Ising nematic-tetragonal transitions. This new scenario has been identified, in particular, 
for 6 = 2, corresponding to stripes of width h = 2 in the ground state. 

In this paper, we perform extensive MC simulations to study the phase transitions in- 



3 



volved in this model for 5 = 2. We present convincing numerical evidence for an Ising 
nematic phase for this coupling and discuss the calculation of the critical exponents in face 
of the constraint related to lattice sizes necessary for observation of the true phenomenology 
of this model. Our study also allows us to figure out how reliable the simulations performed 
with a local update algorithm are. Using histogram reweighting techniques [171 [18] and the 
patching procedure to combine data obtained from simulations at various temperatures, 
we evaluate the specific heat and the order parameter susceptibility over a wide range of 
temperatures, and the free-energy barrier associated with the coexisting phases. This analy- 
sis procedure enables us to explore the existence of distinct maxima as a continuous function 
of the temperature and thus, to estimate the finite-size critical temperatures with high pre- 
cision. Our results show that simulations need to be performed for lattice sizes L at least 
as large as L = 48 so that the Ising nematic phase can be clearly observed. Consequently, 
the conclusions presented in the literature for smaller lattice sizes can be misleading. More- 
over, to ensure the correct determination of the physical quantities and their error bars, one 
needs to have a reasonable number of independent measurements. Thus, to access how the 
long-range dipolar interaction may affect physical estimates, we also evaluate the integrated 
autocorrelation time in energy time series [19]. 

In Sec. 2, we briefly review the histogram reweighting and patching techniques which we 
use to estimate average energies, orientational order parameter and their response functions, 
specific heat and susceptibility, as a function of temperature. Formation of the nematic phase 
is clearly exhibited as we increase the lattice size. In all cases, we estimate the specific heat 
maxima and finite-size critical temperatures by including Jackknife error estimates in the 
analysis |,2Qj. Our histograms for the energy distributions present double-peak structure 
leading to the existence of domain walls between the states of each phase. Thus, we also 
evaluate the free-energy barriers at both phase transitions. The integrated autocorrelation 
time analysis is presented in Sec. 3, which is followed by a brief outlook and our conclusions 
in Sec. 4. 

II. MONTE CARLO SIMULATIONS 

We consider the two-dimensional system described by the Hamiltonian in Eq. Q, where 
the spin variables cr = ±1 are assumed to be aligned perpendicular to the square lattice with 
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N = L"^ sites. We have used the standard Metropohs algorithm to generate configurations 
for lattice sizes L = 16, 32, 48, 56, 64 and L = 72, with 6 = 2. We carried out our simulations 
with periodic boundary conditions to minimize border effects due to the dipolar term in the 
Hamiltonian. Thus, all distances rij must include sites in the infinitely replicated simulation 
box in both directions. 

The boundary conditions add an infinite sum over all images of the simulation box to 
the dipole term of the Hamiltonian. We use Ewald summation technique to compute this 
infinite sum. This technique splits the infinite sum over all images of the system into two 
quickly converging sums: the direct sum, which is evaluated in the real space, the reciprocal 
sum, carried out in the reciprocal space, and a self-interaction correction term |211[22]. We 
have set the Ewald parameter a to 3.5 in all simulations. This parameter determines the 
rate of convergence between the two sums. 

Our MC Markov chains are initialized with random starts. For thermalization, 10^ sweeps 
were discarded for L < 48, and 5 x 10^ sweeps for L = 56, 64 and 72. Our checks have shown 
that after ~ 10^ sweeps the system does not present considerable changes in the energy. A 
detailed analysis about the time evaluation can be found in Ref. [23]. For each temperature, 
measurements rely on 3.4 x 10^ sweeps for L < 64. Our largest lattice size L = 72 relies 
only on 2.7 X 10^ sweeps because even a 10 times larger simulation would not present a 
reasonable number of independent measurements in face of our estimates from the integrated 
autocorrelation time calculations for smaller lattice sizes depicted in Fig. 8. However, the 
patching procedure attenuates this drawback because we simultaneously match simulations 
at different temperatures Tq to obtain the final physical estimates. 

To identify the finite-size critical temperatures, we evaluate the specific heat 



over a continuous range of temperatures by reweighting MC data obtained from simulations 
at Tq. The list of temperatures and our final estimates for the specific heat maxima Cy max 
and their finite-size critical temperatures Tmax is displayed in table 1. All error bar estimates 
rely on 40 Jackknife bins for L < 64 and 20 bins for L = 72. 

With the time series produced at temperatures Tq*'', {i = 1, ■ ■ ■ , P), histogram reweight- 
ing techniques allow us to evaluate physical quantities in neighborhoods of the simulated 




(2) 
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-tO <-^D|max -tmax <-^?;|max -tmax 

16 0.791, 0.830, 0.850, 0.870 3.095(2) 0.8323(2) 

32 0.780, 0.791, 0.812, 0.825, 0.850 3.92(2) 0.7905(3) 

48 0.770, 0.780, 0.791, 0.812, 0.870 3.76(3) 0.7785(7) 4.21(1) 0.8132(2) 

56 0.760, 0.770, 0.780, 0.790, 0.808, 0.860 3.72(4) 0.7727(6) 5.44(2) 0.8084(2) 
64 0.767, 0.773, 0.780, 0.800, 0.808 3.76(5) 0.7726(7) 6.05(3) 0.8045(2) 

72 0.760, 0.770, 0.780, 0.790, 0.807, 0.830 4.05(35) 0.767(2) 5.63(14) 0.800(1) 

TABLE I: Data produced at temperatures Tq and used for patching and reweighting. Tmax^ 
finite-size critical temperatures for the phase transitions observed for system sizes L defined by the 
maximum of specific heat \ max as shown in figures 1 (a) and 1 (b) . 



temperatures, 



1 



A(^)i?. 



(3) 



where nmeas is the total number of measurements in both energy E and physical quantity 
/ time series; A(l/T) = 1/T — I/Tq, and Z{T) is the partition function, 



ZiT) = Y: exp 



A(^)i?. 



(4) 



n=l 

To increase the statistics for the final reweighting, we combine data obtained from inde- 
pendent simulations at temperatures Tq, as shown in table 1, with the condition that the 
P simulations occur at Tq values close to each other [18]. This technique is based on the 



remark that the final estimate /(T), obtained from P independent statistical quantities /j, 
can be calculated as a weighted linear combination 

J{T)=j:w.Jif), (5) 
1=1 

where Wi = Wi{T) are normalized weights, which one takes as the inverse statistical variance 
Wi{T) ^ l/cr^(/j) from each MC time series. The overall constant is determined by the 
normalization condition YliLi ua = ^- This procedure corresponds to increasing the final 
statistics P times. 

Figure 1 shows our final estimates for the averages of specific heat Cy and energy per spin 
(E/N) as a function of temperature. Figure 1(a) displays for L = 16,32, and 48. We 
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FIG. 1: Specific heat Cy and average energy per spin {E/N) as a function of temperature T are 
shown in figures (a) and (b) for L = 16,32,48, and in figures (c) and (d) for L = 56,64,72. The 
inset in figure 1(a) shows Cv with error bar estimates for L = 48. 

have included error bar estimates for L = 48 only in the inset of Fig. 1(a), to have a clear 
presentation of the behavior for different lattice sizes. The numerical simulation with L = 
16 shows a typical striped-tetragonal transition characterized by a peak in Cy. This behavior 
changes dramatically and gives origin to an intermediate thermodynamic transition as we 
increase the lattice size. Moreover, this new transition, the nematic-tetragonal transition, 
presents a more pronounced increase in the specific heat maximum compared with that of 
the striped-nematic transition, which occurs at a lower temperature. 

Figure 2 illustrates the typical spin configurations in each phase for L = 56. This figure 
shows how the directional and translational symmetries are missed as we increase the tem- 
perature toward the tetragonal phase. The corresponding typical spin configurations at the 
transition temperatures are displayed in figure 3. 
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(a) T = 0.710 (b) T = 0.790 (c) T = 0.840 

FIG. 2: Typical spin configurations for L = 56 at temperatures in the: (a) striped, (b) nematic and 
(c) tetragonal phases. The configuration displayed in the striped phase presents E/N = —1.2028, 
Ohv = 0.9873; in the nematic phase, E/N = —1.1172, Ohv = 0.7903; and in the tetragonal liquid 
phase E/N = -1.0201 and = 0.1777. 

The finite-size critical temperatures Tc{L) — T^naxiL) are defined by the maximum of the 
specific heat C^{L). Making Jackknife bins from the combined MC statistics, we reweight 
each bin in order to find the maximum of the specific heat. This procedure leads to final 
estimates of critical temperatures and their error bars in table 1. 

We note that the specific heat peaks occur at lower temperatures as L increases. This is 
another finite-size effect. In fact, finite-size scaling (FSS) arguments predict the following 
behavior for the shift in the critical temperatures 

7;(i^)-r,(oo)-L-v-, (6) 

where v is the critical exponent. This exponent assumes the value 1/d for first-order phase 




(a) r = 0.773 (b) T = 0.808 

FIG. 3: Spin configurations at the transition temperatures for L = 56. The configuration displayed 
in (a) has E/N = -1.1670, = 0.9205; while in (b), E/N = -1.0769 and Ohv = 0.5633. 
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transitions, where d is the system dimension. Another scahng relation for the shift in 
the critical temperatures can be considered in the case of a Kosterlitz-Thouless (KT) type 
transition [21], 

where a and b are constants. These are multiparameter fits, which do not allow us to 
estimate the infinite volume critical temperature Tc(oo) or the proper scenario. In the case 
of a KT type transition, the specific heat peak does not occur at the critical temperature 
but above the transition temperature |25j . 

From the results in table 1 for the finite-size critical temperatures, and because we expect 
to observe both phase transitions also in the thermodynamic limit, it is reasonable to suppose 
that the difference in the critical temperatures will reach a fixed value. 

To further characterize the observed sequence of phases, the degree of orientational order 

mm, 



On. = ^^—^ , (8) 

is also computed for each configuration of the system. The quantities Uh and n„ denote 
the number of horizontal and vertical bonds between oppositely aligned nearest-neighbor 
spins. For the striped ground state, this order parameter takes the value +1, while for 
higher temperatures, in which the orientational symmetry of the striped domains is broken, 
it vanishes. To see how this order parameter and its susceptibility 

x(0,,) = iV((0L)-(0,,)2), (9) 

characterize the striped, nematic, and tetragonal phases, we also study their behavior as a 
function of temperature by means of reweighting techniques. 

We present the maxima of the susceptibilities and their respective finite-size critical tem- 
peratures Tinax in table 2. These temperatures are defined by the maximum of susceptibility 
x{T) as shown in figure 4(a) and (c). Figures 4(b) and (d) illustrate the behavior of the 
order parameter Ohv as a function of temperature for different lattice sizes. They show 
that the break of the directional symmetry happens mainly at the second thermodynamic 
transition temperature Tjp'^ ^ 0.809. The peaks in the susceptibility related to this second 
transition become higher as the lattice size increases. On the other hand, it seems we have 
a trend of lower values for its maximum as the lattice size increases in the case of the first 
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L 




max 


X(^/if)max 


7^(2) 
max 


16 






28.42(2) 


0.8336(2) 


32 






49.4(1) 


0.8258(2) 


48 


17.0(5) 


0.781(2) 


129.7(6) 


0.8140(2) 


56 


13.8(2) 


0.773(1) 


170(2) 


0.8097(2) 


64 


13.5(3) 


0.773(1) 


191(2) 


0.8073(5) 


72 


13.1(9) 


0.768(2) 


209(3) 


0.809(2) 



TABLE II: Susceptibility maxima and the corresponding critical temperatures Tmax as shown in 
figure 4. 

thermodynamic transition temperature TJ^^ ^ 0.77. According to the estimates in table 
2 and their error bars, we are facing an unusual trend of decreasing peaks that continues 
until a limiting value is reached. In any case, this new feature places an extra degree of 
difficulty to the identication of the character of the transition. This situation presents a 
parallel behavior when the maxima of the specific heat is analyzed in table 1. The first 
transition (T^^) ^0.77) seems to present specific heat peaks of roughly the same order of 
magnitude, while the second transition now located at Tjp''^ ^ 0.800 shows a trend of higher 
peaks as the lattice size increases. 

The characterization of the order of the observed thermodynamic phase transitions can 
be achieved by means of FSS analysis. In this context, the specific heat peak is described 
by the relation, 

aUax-iV"/'^^ (10) 

The FSS relation for the peak in the susceptibility is given by 

Xmax~iV^/"^ (11) 

The first-order character is related to the critical exponent dv = 1 [2S1 EZl- From the 
hyperscaling relation a = 2 — and assuming a first-order transition, one obtains ajdv = 1 
and ^jdv = 1. 

Since the phenomenology of this model reveals itself only for rather large lattice sizes, 
our data still present limitations for a reliable FSS analysis. As we are going to show in the 
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FIG. 4: Susceptibility and order parameter for lattice sizes L = 16—72 as a function of temperature. 

next Section, any intensive MC simulations with local update algorithms will suffer severe 
limitations. However, our data allows one to draw some conclusions. 

First of all, let us verify the possible scenario found by the FSS analysis. To this end, 
we must only consider the last four results for the maximum of Cy in table 1 and the 
maximum of xiOhv) in table 2, corresponding to L = 48, 56, 64 and 72. Since the peaks 
for both specific heat and susceptibility related to the first transition seem to converge 
toward a constant value, we may assume a ~ and 7^0. In this case, logarithmic 
corrections should be taken into account to describe their maximum. The supposed values 
for the critical exponents do not support the first-order transition character. If does 
not present a divergent behavior, nor has appreciable finite-size effects, then the hypothesis 
of a Kosterlitz-Thouless type transition must be considered [15|. The hypothesis of a KT 
type transition has been discussed as an alternative to the possibility of the first-order 
character found from MC studies in Ref. |3l [IS]. Here, in spite of the high precision values 
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for the physical estimates, our results for the temperature and response functions may be 
biased when one considers estimates from the two largest lattice sizes, mainly the L = 72 
case, hampering further conclusions. In particular, the lower precision in the specific heat 
at Tj^-* when compared with the results for the second peak at T^^ is related to strong 
autocorrelations in the striped-nematic transition, as presented in the next section. 

The FSS analysis for the second transition yields a = 0.80(2), 7 = 0.77(4) and du = 
1.20(2) for L = 48 — 72. These results clearly support the second order character for the 
nematic-tetragonal phase transition as first suggested theoretically by Abanov [15J, via the 
continuous model, and by Pighm and Cannas [3] from a mean field approach. However, we 
are uneasy about these numerical findings because of their very small godness-of-fit. For 
L = 72, in particular, we may argue that C„ evaluation presents significant bias since its 
value does not follow the main trend, leading the scaling relation to an unlikely fit. However, 
this picture does not change if we discard L = 72. We obtain roughly the same numerical 
estimates for the exponents because Cv{L = 72) has comparatively large error bar. 

Now, we will analyse the energy density distributions and, as explained below, we include 
the calculation of free energies and interface tensions, in our case linear tensions, to gather 
further information about the possible nature of the phase transitions. 

In Fig. 5 we show energy density distributions from data obtained at temperatures Tq. 
Here, we exhibit only energy distributions from simulations at temperatures close to the 
critical ones. In Fig. 5(d) we include two distributions for each lattice size because none of 
the sampled temperatures are close enough to the critical ones to display their true critical 
distributions. Figure 6 shows all histograms of energy distributions for L = 64. 

Let us initially consider the first transition. Figures 5(a) and (c) present histograms 
with clear double-peak structures close to the critical temperatures TJ^^ This is a typical 
first-order phase transition behavior where a latent heat is responsible for the existence 
of those peaks at energies Est{L) and Enem{L), corresponding to the striped and nematic 
phases, respectively. Those states are separated by a minimum at Em{L) corresponding to 
domains describing the coexistence of both phases. Here, we can evaluate the free energy 
barrier AF{L) for this temperature-driven transition. Although it seems clear we face a 
first-order transition, it remains to be seen whether the linear tension is not negligible in 
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FIG. 5: Energy density distributions for lattice sizes L = 48 — 72 from simulations at temperatures 
To. 

the thermodynamic hmit. To this end, we evaluate the free energy [28l 129] 

AF{L) = In ^^(^-*/^'^^) (12) 

at Th, where the finite lattice critical point T = Th{L) is defined such that the histogram 
peaks have equal heights Pii^Est/N, Th) = PL{Enem/N, Th). We proceed performing patching 
and reweighting of histograms [30], in order to evaluate AF{L) at both transitions. This 
leads to the results in Fig. 7 for both transitions. In this figure, the solid and dashed lines 
are the linear extrapolation in 1/L, 

where /o is the linear tension and /i stands for possible FSS correction. Here, we assume the 
two-gaussian approximation [281 E] to describe the coexisting phases in this model. This 
procedure yields /o = 0.019(4) and /i = 0.37(16) for the striped-nematic transition. 
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FIG. 6: Energy density distributions for L = 64 obtained from simulations at temperatures Tq. 

Now, looking at the normalized energy density distributions for the second transition, 
Fig. 5(b) and (d), we also observe the characteristic double-peak behavior although the 
histograms in Fig. 5(b) have less pronounced peaks. We observe higher peaks for L = 64 
and 72 after reweighting to the critical point where PiiEnem/N, T^) = PL^Etetra/N, T^). For 
this transition one obtains /o = 0.0077(2) and /i = —0.269(23). The second transition is 
clearly weaker compared with the first one. Just to place some comparative figures to clarify 
the meaning of these numbers, one has /o = 0.04735 for 2D Potts model with 10 states and 
/o = 0.010395 for this model with 7 states. These numbers are exact results. The first-order 
character of the model with 10 states is easily identified and is sufficient to work with lattice 
sizes up to 50. It is well known that this size is too small for such identification in the 
case of g = 7. 

Long-range interactions always impose strong finite-size effects. Although our analysis 
procedure based on reweighting techniques with multiple histograms extracts the maximum 
information available in the MC data, some aspects still prevent us from definitive conclu- 
sions about the nature of the phase transitions. In all cases, the bulk correlation length 
^ plays an important role. Actually, ^ and the extension of finite systems L are the only 
relevant length scales. The validity of the FSS ansatz requires that L > C,- The first tran- 
sition seems to be of KT type or even of first-order nature. For a KT type transition, it is 
predicted that ^ has the following behavior [32] at Tc 





^ ~ exp(6/t 




(14) 
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FIG. 7: Extrapolation of the free energy as a function of L^^ for the striped- nematic ( — ) and 
nematic-tetragonal ( ) transitions. 

where t = T/Tc — 1 is the reduced temperature. For first-order transitions ^ remains fi- 
nite and for second-order ones it behaves as ^ ~ t^". Thus, at Tc the correlation length 
diverges exponentially, while divergence is the power law in a second-order phase transition. 
Unfortunately, a dedicated work would be required to clarify such scenario [33j. If present, 
the KT type transition would be responsible for eliminating the long-range order in the 
striped phase through the formation of unbound dislocations [15] • The reason for consid- 
ering a KT type transition is the absence of divergences or significant finite-size effects in 
the response functions. However, this absence may be justified by the enormous difficulty 
in sampling configurations at this transition, as revealed later in our analysis of integrated 
autocorrelation times. 

Our results from the interface tension analysis for the nematic-tetragonal transition indi- 
cate a very weak first-order phase transition. This transition has been described as second- 
order from the FSS analysis of the response functions specific heat and susceptibility. There- 
fore, the first-order character exhibited by the histograms for the larger lattice sizes L = 64 
and 72 does not seem to be caught by this FSS analysis. This implies that the fiuctuations 
in the energy and order parameter are somehow damped in our data. This may result from 
the small lattice sizes we work with, where the true first-order character is not effective and, 
presumably, combined with the need of a larger number of independent configurations in 
both phases. Actually, according to our results the effect of the free-energy barrier should 
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be more severe between the striped and nematic phases. In turn, rehable calculation of the 
free-energy barriers requires enough configurations at Em that are exponentially suppressed 
by the Boltzmann factor in canonical MC simulations. 

III. INTEGRATED AUTOCORRELATION TIME 

Local MC updates are not efficient, and the produced data correspond to successive states 
of a Markov chain. This chain may be highly correlated, introducing bias to the estimates 
of the physical quantities and unreliable variances if the sampling is not large enough. Also, 
related to systems with long-range interactions, the presence of a very slow approach to 
equilibrium introduces a further degree of complexity. This has been highlighted in a recent 
study by Cannas et al. [23j on the dipolar model. They found stronger metastabilities 
associated with the striped-nematic transition compared with the nematic-tetragonal one. 
Here, we complement our study by evaluating the integrated autocorrelation time, to assert 
the effectiveness of the collected data at each temperature Tq with local update algorithms. 
To this end, we define the normalized autocorrelation function for the energy time series 

a\E) 

where (t{E) is the energy standard deviation. The integrated autocorrelation time Tint, 

1 rir 

r^nt = .P(0) + E P(0 (16) 

estimates the number of independent data points in a long sequence of rir MC measurements. 
Actually, twice this value is the proper quantity to effectively evaluate this independence. It 
is usually expected that the autocorrelation function decays as a power law. However, this 
behavior is changed by a temperature dependent factor near critical points [12] . 

Our calculations of 2Tint are shown in Fig. 8 for lattice sizes L = 16, 32, 48 and 56, where 
we have introduced the notation k to describe the sample lengths, Ur = 2^. 

We obtain 2rj„t = 800 measurements (L = 16) for temperatures close to the finite-size 
critical point = 0.832 and smaller values for temperatures in the tetragonal phase. For 
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FIG. 8: Integrated autocorrelation time 2Tint for lattice sizes L = 16 (a), L = 32 (b), L = 48 (c) 
and L = 56 (d) as a function of time series length 2'^. 

the L — 32 case, Fig. 8(b) shows a larger value: 2Tj„t = 6 x 10^ for T — 0.780 in the striped 
phase and 4 x 10^ for T — 0.791, in the vicinity of the specific heat maximum. Again, as 
expected, smaller values are observed for samples in the tetragonal phase. 

Figure 8(c) displays the results for L = 48, whose specific heat data present two critical 
temperatures. For the temperature inside the striped phase (T = 0.770), one obtains a 
larger integrated autocorrelation time compared with the one achieved at the striped-nematic 
transition 2Tint{T — 0.780) = 5 x 10^. In the nematic phase, one obtains 2Tint{T — 0.791) = 
4 X 10^. By increasing the temperature toward the next critical point (Tc = 0.8132), one 
obtains 2rj„f(T = 0.812) = 2.4 x 10^. As the temperature increases further, it is observed a 
strong decrease in 2rint, reaching values as small as 45 for T = 0.870. 

Figure 8(d) shows huge numbers for this autocorrelation time, larger than 2 x 10® for data 
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collected in the striped and at the striped-Ising nematic transition. In the Ising nematic 
phase, one has 2Tint{T = 0.790) = 9.6 x 10^. At the Ising nematic-tetragonal transition one 
obtains 2Tint(T = 0.808) = 5 x 10^, and a smaller value 60 at the temperature T = 0.860 is 
achieved in the tetragonal phase. 

Here, we could try to describe autocorrelation times for the energy observable as a simple 
power law 2Tint,E = AL^"'\ where Zint is the associated dynamic critical exponent. A rough 
estimate gives Zm ~ 6.2 and Zint ~ 4.6 for the striped-nematic and nematic-tetragonal phase 
transitions, respectively. These are large numbers since local MC algorithms generally have 
a dynamic critical exponent z ^ 2 [34J. For non-local update algorithms, such as the 
Swendsen-Wang cluster algorithm, this exponent is further reduced to z = 0.222(7) in the 
pure 2D Ising model. 

From the autocorrelation time analysis, we may infer the existence of very long-lived 
states near and at the striped-nematic transition. These states may lead to stronger bias in 
the determination of the response functions at this transition than at the nematic-tetragonal 
transition, where Tint are orders of magnitude smaller. 

IV. CONCLUSIONS 

Our systematic analysis by reweighting techniques in multiple histograms describes the 
existence of two phase transitions, which correspond to a two-step melting process leading to 
a disordered state. This long-range interaction model is described by a complex phase dia- 
gram with lines corresponding to striped-nematic and nematic-tetragonal phase transitions, 
depending on the coupling 6. We have seen that the true phenomenology of this model can 
be observed for large lattice sizes only, and that the controversial nematic phase is detected 
in a narrow range of temperatures. As 6 increases, one observes larger stripe width. As 
realized for h = 2, we also expect that for larger values of h, it will be necessary to increase 
further the system size to obtain reliable results. This can be explained by the need of 
inserting dislocations in the striped-domain structure. The formation of bound dislocation 
pairs proliferates as the temperature increases, leading to a phase transition. 

Our numerical results clearly show the locations of the phase transitions, which are 
strongly lattice size dependent. Moreover, the response functions also have strong finite- 
size effects, which frustrate any convincing evaluation of the critical exponents from simple 
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FSS analysis. 

The histogram analysis has shown a strong first-order character for the striped-nematic 
transition. However, this character may introduce strong metastabilities, leading to bias in 
our simulation. This is because the sampling may spend most of the time in one of two 
phases, mainly in the case of our largest lattice sizes. According to our integrated autocor- 
relation time calculations, we have learnt that to obtain 10'^ independent configurations for 
L = 72, about 10^*^ and 10^ MC sweeps are required for the striped-nematic and nematic- 
tetragonal transitions, respectively. We have attenuated this limited sampling by matching 
simulations at different temperatures but close to each other. However, it is not just a matter 
of how long one needs to perform the simulations. Our analysis has also revealed the need 
for larger lattice sizes, so that the character of the phase transitions is exposed. It is likely 
that the lattice sizes are not large enough compared to the correlation lengths. This because 
our critical exponent estimates from FSS indicate continuous transitions, in contrast with 
the results for the free-energy barriers. The asymptotic behaviors of the free-energy barriers 
between domains indicate that the interface tensions do not vanish in the thermodynamic 
limit. 

Another controversial point still remains, the possibility of a KT type transition for the 
nematic-tetragonal transition. Our results for the response functions do not discard this 
possibility. 

The integrated autocorrelation time calculations show stronger critical slowing down in 
both phase transitions when compared to the pure 2D Ising model. From this analysis we 
may conclude that there are very long-lived states near and at the striped-nematic phase 
transition as opposed to the nematic-tetragonal transition, where Tint present smaller orders 
of magnitude. These results are very important to figure out how reliable the MC sampling 
for this system with long-range interaction is. Moreover, the use of important techniques for 
MC data analysis establishes the limits one can reach with local update algorithms. More 
efficient results to obtain a deeper understanding of the complex striped-nematic-tetragonal 
phase transitions seem possible with a variant of WoUf cluster algorithm [35] or by means 
of generalized ensemble methods. 
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